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We investigate the critical dynamics of the Hybrid Monte Carlo algorithm approaching the chiral limit of 
standard Wilson fermions. Our observations are based on time series of lengths 0(5000) for a variety of observ- 
ables. The lattice sizes are 16^ x 32 and 24^ x 40. We work at /3 = 5.6, and k = 0.156, 0.157, 0.1575, 0.158, with 
0.83 > > 0.55. We find surprisingly small integrated autocorrelation times for local and extended observables. 
The dynamical critical exponent z of the exponential autocorrelation time is compatible with 2. We estimate the 
total computational effort to scale between and 1/^4 towards the chiral limit. 



1. INTRODUCTION 

Considerable effort has been spent since the 
early days of exact full QCD simulations with Hy- 
brid Monte Carlo methods [Q to optimize both 
performance and de-correlation efficiency. The 
main issues are parameter tuning, algorithmic im- 
provements and assessment of de- correlation. 

HMC parameter tuning has been the focus of 
the early papers. At this time one was restricted 
to quite small lattices of size w l** . As main 
outcome, trajectory lengths of 0(1) and accep- 
tance rates of « 70 % have been recommended. 

The last three years have seen a series of im- 
provements: Inversion times could be reduced by 
use of the BiCGstab algorithm and parallel 
SSOR preconditioning [Q. Together with chrono- 
logical inversion M, gain factors between 4 and 8 
could be realized 

Less research has been done on the critical dy- 
namics of HMC. In Ref. [Q, the computer time 
required to de-correlate a staggered fermion lat- 
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tice has been estimated to grow like T oc niq* . 
A similar law has been guessed in Ref. (ll| with 

-23 

T oc m^^ for two flavours of Wilson fermions. 
However, so far we lack knowledge about the scal- 



ing of the critical behaviour of HMC: a clean de- 
termination of the autocorrelation times requires 
by far more trajectories than feasible at the time. 
Ref. quotes results in SU(2), the APE group 
found autocorrelation times in the range of 50 
for meson masses [|l^, the SCRI group gives esti- 
mates for Tint between 9 and 65. However, they 
seem to behave inconsistently with varying quark 
mass [p^ . 

SESAM and T^L have boosted the trajectory 
samples to 0(5000), i.e., by nearly one order of 
magnitude compared to previous studies . The 
contiguous trajectories are generated under stable 
conditions and allow for reliable determinations of 
autocorrelation times from HMC simulations. 

2. AUTOCORRELATION TIMES 

The finite time-series approximation to the true 
autocorrelation function for At, t ^ I, . . . , Imc is 

tMC—t /tMC \ 

E A^A^^t-j^AEAA 

C'^{t) = ^ ^^^i (1) 

tMC — t 

The exponential autocorrelation time is the in- 
verse decay rate of the slowest contributing mode 
with p'^it) being C^{t) normalized to p"^(0) = 1, 
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Figure 1. p{t) at Ksca — 0.157 for Dixi, and A. 



the length of the thermalization phase of the 
Markov process. To achieve ergodicity the sim- 
ulation has to safely exceed sup^ {r^p}. The in- 
tegrated autocorrelation time reads: 



tMC- 



(2) 



In equilibrium r^^j characterizes the true statis- 
tical error of the observable A. The variance 
a\ = 2Tj^(CTg is increased by the factor 2t^^ com- 
pared to the result over a sample of N indepen- 
dent configurations, r^^j is observable dependent. 

3. RESULTS 

The results given here can only be a small ex- 
cerpt of our investigation to be presented else- 
where in more detail p^ . The autocorrela- 
tion is determined from the plaquette (Dixi) 
and extended quantities like smeared light meson 
masses, ttItt and nip, and smeared spatial Wilson 
loops (Daxs and Digxie)- They exhibit a large 
ground state overlap per construction. Further- 
more we exploited the inverse of the average num- 
ber of iterations A = N^^y of the Krylov solver 
which is related to the square root of the ratio 
of the minimal to the maximal eigenvalue of the 
fermion matrix 

We illustrate the observable dependency of the 
autocorrelation in Fig. |l|. Using 0(3000) trajecto- 
ries we achieve a clear signal for the autocorrela- 
tion function. This is the case also for the dynam- 
ical samples at Kgca = 0.156, 0.157 and 0.1575 as 
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Figure 2. Kgoa dependency of A on the 16"^ x 32 
lattice. 



well as on the 24r x 40 lattice with Ks, 



0.1575 



and 0.158. A appears to give an upper estimate to 
the autocorrelation times of 'fermionic' and 'glu- 
onic' observables measured^. 

Shifting the sea-quark mass towards the chiral 
limit, we observe increasing Texp and Tint- Fig. ^ 
sketches this dependence of on the 16^ x 32 
lattice. 

The volume dependence of the autocorrelation 
was compared for two different lattice sizes at 
equal Ksca = 0.1575. Unexpectedly, we found Texp 
and Tint, measured in units of trajectory num- 
bers, to decrease by about 50 % for Dixi and 
30 % for A, while switchng from the 16"^ to the 
24^^ system. As we chose the length of the HMC 
molecular dynamics as T = 1 on the 16'^ x 32 sys- 
tem and T = 0.5 on the 24^ x 40 lattice, these 
numbers are even more surprising. The origin of 
this phenomenon needs further investigations. 

A compilation of Texp is given in Fig. the 
values for Tint being similar. 

The quality of our data allows to address the 
issue of critical slowing down for HMC. The ap- 
proach to Kgca amounts to a growing pion cor- 
relation length = l/rriT^a. The autocorrela- 
tion time is expected to scale with a power of 
S,, T — e^^. The dynamical critical exponent z 
governs the scaling of the compute effort, while 
the knowledge of e allows to assess its absolute 
value. Our results are given in Tab. |l|. For the 



^Slower modes might exist for the topological charge 



3 



Table 1 



Critical exponents z. The 16^ x 32 results are based on full, the 24^ x 40 numbers on partial statistics. 



Size 


16^ X 32 


24^ X 40 


Observable 


□ixl 
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^8x8 


□ixl 


A 


1^16x16 




1.8(4) 


1.5(4) 


1.2(5) 


2.2 
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0.1 




1.4(7) 


1.3(3) 


1.3(3) 
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Figure 3. Exponential autocorrelation times. 



smaller lattice they tend to be below z = 2, albeit 
they are still compatible with 2. It seems that 
extended observables on the large lattice yield 
smaller values. However, we have to wait for 
larger samples to arrive at conclusive answers. 

Finally we try to give a conservative guess for 
the scaling of the compute time T required for 
de-correlation. The maximal length of has 
been limited to V^/^tt « 4 in order to avoid fi- 
nite size effects. With fixed, the volume factor 
goes as Furthermore we found A oc ^^^'^ for 
BiCGstab. z lies between 1.3 and 1.8, taking the 
result from the 16'^ system. 

In order to keep the acceptance rate constant, 
we reduced the time step from 0.01 to 0.004 with 
increasing lattice size, while we have increased the 
number of HMC time steps from 100 to 125. 

At the same time, as we have already men- 
tioned, the autocorrelation time for the worst case 
observable A went down by 30 % compensating 
the increase in acceptance rate cost! In a conser- 
vative estimate, we thus would assume that the 
total time T scales as to This guess trans- 
lates to T oc ^2 - y2.2_ 



4. SUMMARY 

The autocorrelation times from HMC under 
realistic conditions are smaller than anticipated 
previously, staying below the value of 60 trajec- 
tories for Texp- z is compatible with 2. The 
computer time to generate decorrelated configu- 
rations scales according to , a very encouraging 
result compared to previous estimates . 
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